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NUMERICAL SOLUTION OF THE TIME-FRACTIONAL 
FOKKER PLANCK EQUATION WITH GENERAL FORCING* 

KIM NGAN LEt, WILLIAM MCLEANt, AND KASSEM MUSTAPHA* 

Abstract. We study two schemes for a time-fractional Fokker-Planck equation with space- 
and time-dependent forcing in one space dimension. The first scheme is continuous in time and is 
discretized in space using a piecewise-linear Galerkin finite element method. The second is continuous 
in space and employs a time-stepping procedure similar to the classical implicit Euler method. We 
show that the space discretization is second-order accurate in the spatial L 2 -norm, uniformly in time, 
whereas the corresponding error for the time-stepping scheme is 0(fc“) for a uniform time step k, 
where a G (1/2, 1) is the fractional diffusion parameter. In numerical experiments using a combined, 
fully-discrete method, we observe convergence behaviour consistent with these results. 

Key words. Time-dependent forcing, finite elements, fractional diffusion, stability, Gronwall 
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1. Introduction. We investigate the numerical solution of the inhomogeneous, 
time-fractional Fokker-Planck equation [10], 

Ut - Kadl~°^Ua:x + Ma ^ = 9, (1-1) 

for 0 < X < L and 0 < t < T, with initial data m(x, 0) = uo(x) and subject to 
homogeneous Dirichlet boundary conditions u{0,t) = 0 = u{L,t). (We use subscripts 
to indicate partial derivatives of integer order with respect to x or t; for instance, 
Ut = du/dt.) The parameter Ua is the generalized diffusivity constant, is the 

generalized friction constant, and the driving force F and the source term g are 
permitted to be functions of both x and t. The subdiffusion parameter a satishes 
0 < a < 1, and the fractional time derivative is interpreted in the Riemann-Liouville 
sense; thus, where /“ is the fractional integral of order a, 

pt 

I°‘v{t) = / OJait — s)v{s) ds with OJctit) = . . . 

Jo r(Q;) 

In 1999, Metzler et al. [15] used a discrete master equation to model the behaviour 
of subdiffusive particles in the presence of a driving force F(x), showing that in the 
diffusive limit the probability density u(x, t) for a particle to be at position x at time t 
obeys a fractional Fokker-Planck equation of the form 

Ut - dl~°‘{KaUxx - ^J.a^{Fu)J =0. (1.2) 

Subsequently, Henry et al. [10] considered the more general case when F = F{x,t) 
may depend on t as well as x, and showed that u obeys (1.1) with g = 0. The two 
equations coincide if F is independent of t, but if the forcing is time-dependent then 
(1.2) does not properly correspond to any physical stochastic process [9]. 
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Various numerical (time stepping finite difference) methods have been proposed 
for solving (1.2), usually for F assumed to be either a constant or a function of x only. 
The starting point was often to rewrite equation (1.2) in the form 

^l~^{Fu)^=Q, (1.3) 

in which the first term is a Caputo fractional derivative. Indeed, (1.3) is in some ways 
more convenient than (1.2) for constructing and analyzing the accuracy of numerical 
schemes. However, the simpler form (1.3) is not applicable in our case because F may 
depend on t. 

For the numerical solution of (1.3) with F = F{x), Deng [6] transformed the 
equation to a system of fractional ODEs by discretizing the spatial derivatives and 
using the properties of Riemann-Liouville and Caputo fractional derivatives, and then 
applied a predictor-corrector approach combined with the method of lines. This work 
also presented a stability and convergence analysis. Cao et al. [3] adopted a similar 
approach for (1.2) and solved the resulting system of fractional ODEs using a second 
order, backward Euler scheme. Chen et al. [4] studied the stability and convergence 
properties of three implicit finite difference techniques, in each of which the diffusion 
term was approximated by the standard second order difference approximation at the 
advanced time level. In related work, Jiang [11] established monotonicity properties 
of the numerical solutions obtained by using these schemes, and so showed that the 
time-stepping preserves non-negativity of the solution. Based on this property, a new 
proof of stability and convergence was provided. 

Fairweather et al. [8] investigated the stability and convergence of an orthogonal 
spline collocation method in space combined with the backward Euler method in 
time, based on the LI approximation of the fractional derivative. In an earlier work, 
Saadmandi [17] studied a collocation method based on shifted Legendre polynomials 
in time and Sine functions in space. Recently, Vong and Wang [18] have analysed a 
high order, compact difference scheme for (1.3), and Cui [5] has considered a more 
general fractional convection-diffusion equation, 

- {aux)x + bux + cu = f, 

with coefficients a, 6, c that may depend on x and t, applying a high-order approxi¬ 
mation for the time fractional derivative combined with a compact exponential finite 
difference scheme for approximating the convection and diffusion terms. Stability 
(using Fourier methods) and an estimate for the local truncation error were obtained 
in the case of constant coefficients. We are not aware of any previous analysis on the 
numerical solution of (1.1) for a general F depending on both x and t. 

In Section 2 we gather together some preliminary results needed in our subse¬ 
quent analysis, including continuous and discrete versions of a generalized Gronwall 
inequality involving the Mittag-Leffler function in place of the usual exponential. One 
of these results (Lemma 2.4) holds only for 1/2 < a < 1 so much of our theory requires 
this restriction. Section 3 deals with a spatial discretization of (1.1) by a continuous, 
piecewise-linear Galerkin finite element method. We prove stability of the scheme in 
Theorem 3.2 and, under weaker assumptions on F but with a worse bound, in The¬ 
orem 3.3. An error estimate follows in Theorem 3.4 showing second-order accuracy 
in Loo(( 0,T),L2(D)), where fl = (0, L) denotes the spatial interval. We then study 
a time stepping scheme in Section 4, proving a stability estimate in Theorem 4.3 and 
then an error bound in Theorem 4.4, assuming a constant time step k. This scheme, 
which is continuous in space, is formally first-order accurate but, owing to the weakly 
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singular kernel in the fractional integral, we are able to show only that the error in 
L 2 (n) at the nth time level is 0(fc“). Section 5 reports on numerical experiments with 
a fully discrete scheme based on the semi-discrete ones analyzed in Sections 3 and 4. 
We observe 0{k + h^) convergence when a is close to 1, or when we use an appropri¬ 
ately graded mesh in time. The experiments give no evidence that the methods fail 
ifO < a < 1/2, although the convergence rate deteriorates as a decreases when using 
a uniform time step. We also apply our method to a problem from a recent paper by 
Angstmann et al. [1] and investigate whether the regularity of the initial data affects 
the stability of the methods. A brief appendix proves a technical result (Lemma A.2) 
used in showing stability of the time-stepping procedure. 

2. Technical preliminaries. Lemmas 2.1-2.4 below summarize some proper¬ 
ties of fractional integrals that will be needed in our analysis. In each case, we assume 
that the function v{t) is defined for 0 < t < T and takes values in a Hilbert space 
with inner product (•, •) and norm || ■ ||, and is sufficiently regular for the integrand 
on the right-hand side to be absolutely integrable. 

Lemma 2.1. For 0 < a < 1 and 0 < t < T, 


1 - a Jo 


Proof. Put w(t) = so that v{t) — 'u(O) = I^Vt = and 

||-i;(t)-u( 0 )f < wi_a/ 2 (t - 5)11^(5)11 ds 


< 


giving the desired bound, because r(l — a/2) > 1. □ 

Lemma 2.2. For 0 < a < 1, 

/ dt < - ^ j {rv{t),v{t))dt. 

Jo cos(a7r/2) Jq 

Proof. Mustapha and Schotzau [16, Lemma 3.1 (ii)]. □ 

Lemma 2.3. For 0 < a < 1, 


\\I°‘v{t)\\'^ dt < u}a+i{T) [ uJa{T-t) [ \\v{s)\\'^ dsdt. 

Jo Jo 


Proof. Since 


ii/“^^(i)ir < 


ft 2 

LOait - s)||'y(s)|| ds 


’0 


- (^J ^ait-s)ds^ Wa(t - s)||?;(s)|p d5^ 
= U}a+l{t) Wa(t - 5)||u(5)|pds 

Jo 
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we have 


cT pt 


\\I^v{t)fdt<uj^+iiT) 


uJa{t — s)||t;(s)||^ ds dt, 


'0 JO 


and the double integral on the right equals 


[ lk(®)lP [ 0 Ja{t — s) dt ds = [ ||?^(s)|p f UJa{T — t) dtds. 
Jo Js Jo j S 

The result follows after reversing the order of integration again. □ 

Lemma 2.4. For 1/2 < a < 1, 

Proof. The identity dl~°‘v{t) = v{0)u!a(t) + I°‘vt(t) implies that 

r \\dl-‘^vm^dt<2\\v{0)\\^ ru;^{t)^dt + 2 \\rvt\\^ dt, 

Jo Jo Jo 

and the Cauchy-Schwarz inequality gives 

T T / t \ 2 

j \\I°'vt\\‘^dt< J a;a(< - s)||ut(s)|| dsj dt 

J yj U}a{t - sf dsj yj ||ut(s)||^ dsj dt, 


< 


so it suffices to note that f* uJa(t — s)^ ds < ((2a — l)r(a)^) for a > 1/2. □ 

The existence and uniqueness of our spatially discrete solution to (1.1) will follow 
from the following result for an m x m system of weakly singular integral equations. 
Here, | ■ | may denote any matrix norm on induced by a norm on R™. 

Theorem 2.5. There exists a unique continuous solution u : [0, oo) R™ to the 
linear Volterra integral equation 

u{t) + f K{t, s)u{s) ds = g{t) for0<t<oo, 

Jo 

if the following conditions are satisfied: 

1. g : [0, oo) R™ is continuous; 

2. K[t, s) € R™^™ is continuous for 0 < s < t < oo; 

3. for any continuous function v : [0, oo) R'”, the integrals 


f K{t,s)v(s) ds and j \K{t,s)\ds 
Jo Jo 


exist and are continuous for 0 < t < oo; 

4- there exist constants 7 > 0 and e > 0 such that 


f e ^^\K{t,s)\ds < 1 — e for0<t<oo. 
Jo 
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Proof. Becker [2, Corollary 2.3]. □ 

Our stability analysis of the spatially discrete solution makes use of the following 
weakly singular Gronwall inequality, involving the Mittag-Leffler function 


Ep{z) = 


OO 

^ rXT 

n=0 ^ 


r/3)’ 


( 2 . 1 ) 


The usual Gronwall inequality is just the special case /3 = 1, because Ei{z) = e^. 

Lemma 2.6. Let /3 > 0 and T > 0. Assume that a and b are non-negative and 
non-decreasing functions on the interval [0,T]. If y : [0,T] —> R is a locally integrable 
function satisfying 

0 < y{t) < a(t) + b(t) f ujp{t — s)y{s) ds for 0 <t < T, 

JQ 

then 

y{t) < a{t)Ep{b(t)t^^ for 0 <t <T. 


Proof. Dixon and McKee [7, Theorem 3.1]; Ye, Gao and Ding [19, Gorollary 2]. □ 
We also use a discrete version of this Gronwall inequality to establish stability of 
our time stepping procedure. 

Lemma 2.7. Let 0</3<l, 7V>0, fc>0 and tn = nk for 0 < n < N. Assume 
that {An)n=o ** ® non-negative and non-decreasing seguence, and that B > 0. If the 
sequence {y^)n=o satisfies 

n—l 

0 < y" < + Bk ^ ujp{tn — tj)y^ for 0 < n < N, 

j=o 


then 


2 /” < AnEp{Bt^) for0<n< N. 


Proof. Dixon and McKee [7, Theorem 6.1]. □ 

3. Spatial discretization. We choose a partition 0 = Iq < 2 ;i < X 2 < • • • < 
xp = L oi the spatial interval II = (0, L) and denote the length of the pth subinterval 
by hp = Xp — Xp-i for 1 < p < P. With h = maxi<p<p hp, we define the usual 
space S/i of continuous, piecewise-linear functions that satisfy the Dirichlet boundary 
conditions, so that C 77^(0). Recall that F = F{x,t) and g = g{x,t). In our 
notation, we will often suppress the dependence on x and think of u = u{x, t) as a 
function of t taking values in L 2 {Il). We also assume that «„ = /io, = 1. 

In the usual weak formulation of (1-1), we seek u satisfying 

{ut,v)-I {dl~°‘uxjVx) - {Fdl~°‘u,Vx) = {g{t),v) for u € i7o(f^), (3.1) 

with u(0) = Uq, where (•, •) denotes the inner product in L 2 {Il). For our error analysis, 
it is useful to consider a slightly more general, spatially discrete version of (3.1), in 
which Uh : [ 0 ,r] satisfies 

{uht,v) + {dl~°‘uhx,Vx) - {Fdl~°‘uh,Vx) = {g{t),v) + {g*{t),Vx) (3.2) 
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for V G Sh, with Uh{0) = u^h where Ug « Ugh £ ^h, and where Uht = duh/dt. Thus, 
if gt{x,t) = 0, then Uh is the standard Galerkin finite element solution of (3.1). 

To show the existence and uniqueness of satisfying (3.2), define the linear 
operator Bh{t) : §/i —>• S/j (which depends on t through F) by 

{Bh(t)v,w) = {vx,Wx) - {Fv,Wx) for v, w G Sh, 
and the finite element function gh{t) £ §/t by 

{gh{t),w) = {git),w) + {g*{t),Wx) for w £ S/^. 

The variational equation (3.2) is then equivalent to 

Uht + Bh{t)dl~°‘uh = gh{t). 


Integrating with respect to t, we find that Uh satisfies the Volterra equation 
'^h{t) + [ Kh{t, s)uh{s) ds = Gh{t) for 0 < t < T, 
with the weakly-singular kernel 

Kh(t, s) = Bh{t)u)a{t -s)- Bht{T)uJa{T - s) dr 

J S 


and right-hand side 

Ghit) =ugh+ / gh{s) ds. 
Jo 


Theorem 3.1. If F € W))q((0,T);Too(^^)) and g, g^, G Li(^(0,T); L 2 (fl)) ^ then 
for any ugh G Sh there exists a unique continuous Uh '■ [0, oo) —>• satisfying (3.2) 
for all V GSh, with Uh(0) = Ugh- 

Proof. Let | • | denote any norm on the finite dimensional space Sh- Our assump¬ 
tions on F, g and ensure that Gh satisfies condition 1 of Theorem 2.5, and that 
Kh satisfies conditions 2 and 3 (after fixing any basis for Sh). Furthermore, 


\Kh{t,s)\ < CF,h(^UJa{t - s) + UJair - s) dx^ = CF,h[aJa{t - s) +UJi+a{t - s)] , 


and, denoting the Laplace transform by £, 


SO 


pCio 

/ — s) ds < / e~^^uJa{s) ds = Cujai'y) = "f~ 

Jo Jo 

f s)| ds < + T”'”"], 

^0 


and condition 4 follows for 7 sufficiently large. □ 

Theorem 3.1 gives no meaningful stability result for Uh{t) (because Gf.h from the 
proof grows rapidly as h — 0 ) but an energy argument yields the following estimate. 
We use the abbreviation ||?;||r for the norm in 

Theorem 3.2. If, in addition to the assumptions of Theorem 3.1, 
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1. F^{x, t) > 0 for 0 < X < L and 0 < t < T; 

2. 1 + F{x, tY < Cp for 0 < X < L and 0 < t < T; 

3. 112 <a< 1; 

then 

\\uh{t) - UohW^ < ( 3 ^^ 115 ( 5 ) 1 ^ + ||5*(s)lP) ds 

+ """'■ii- /■>''»< 

Proof. Using (3.2), we find that the function Wh = Uh — uo/i satisfies 

(wht.v) + {dl~°'whx,Vx) - {Fdl~°‘wh,Vx) = { 9 {t),v) + {J{t),Vx) 

for all V G S/i, where J{x,t) = g^{x,t) + F{x,t)dl~‘^uoh{x) — {uoh)xix,t). Choos¬ 
ing V = dl~°‘wh(t) G Sh, 

{wht,dl~°‘wh) + \\dl~°‘whx\\'^ - {Fdl~°^Wh,dl~°^Whx) 

= {g{t),dl-<^Wh) + {J{tldl-<^wux), (3.3) 

and since dl~°‘wh{x,t) = dl~°‘uh{x^t) — ijJa{t)uQh{x) = 0 if a; G {0,L}, integration by 
parts gives 

{Fdl-^Wh.dl-‘^Whx)= ! F\{{dl-‘^Whf)jx = - j Fx\{dl-‘^Whf dx. 

Jo Jo 

Hence, by assumption 1, 

{wht,dl~°‘wh) + \\dl~°‘whx\\'^ < {g{t),dl~°‘wh) + {J{t),dl~°‘whx), (3.4) 

and the Poincare inequality, ||r!|p < for v G i7o(U), implies that 

{g{t),dl--wn) < ^\\9{t)r + l\\dl-^w^,xr. 

Using {J{t),dl~°‘whx) < 3ll>/(0lP + \\\dl~°'whx\\^, and noting that dl~°‘wh = F^wm 
because WhiO) = 0, it follows from (3.4) that 

{wht,I°'wht) = {wht,dl-°‘wh} < ^||5(0lP + 

By Lemmas 2.1 and 2.2, and using the inequality cos(Q;7r/2) >1 — 0 , 

\\wh{T)\\'^<-^ -^ / {l°‘wht,Wht)dt 

{l-ayjo 

rpl—a /r2 rT t \ 


and since [d] °^UQh){x,t) = uJait)uQh{x), the Cauchy-Schwarz inequality gives 

lk(f)|| < ||5*ll + \/^u;a(<)||Mo/i||i- 
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Hence, by assumption 2, 

i / ||J(<)f (it < / Wg^tWdt + CFWuohWl [ ujaitfdt, (3.5) 

^ Jo Jo Jo 

and assumption 3 means that uJa(t)^ dt < dt = j(2a — 1), implying 

the desired estimate. □ 

For applications, the condition > 0 seems unnaturally restrictive. In the next 
result, we show that it is not necessary for stability, but the resulting bound grows 
more rapidly with t, owing to the use of the weakly singular Gronwall inequality. 

Theorem 3.3. If we drop assumption 1 from the hypotheses of Theorem 3.2, 
then 

\\uh{t) - Uohf < ^ ^ (i||c,(s)f + ||g*(s)f) ds 

+ 2 J- 1 ll^o(i|li^ forO<t<T. 

Proof Recall that dl~°‘wh = I°'Wht because Wh{0) = 0, so (3.3) implies that 
{rwuuwut) <\\\FI^WMf + k\\9(t)f + + \\\J{t)f. 

By Lemma 2.2, 

yh{T)= [ dt < -^— [ (l°‘wht,Wht) dt, 

Jo 1 - a Jo 

so if we let 

T 

a(T) = 2(Y^^ iimr + wmndt 

then 

yh{T) < a{T) + r Wl^wutf dt. 

2(1 — a) Jq 

Since I°‘wht = Lemma 2.3 implies that 

yhiT) < a{T) + b{T) [ uj^/ 2 {T - t)yh{t) dt where b{T) = 

Jo 2(1 - a) 

Hence, using Lemma 2.1 and the Gronwall inequality of Lemma 2.6, 

lk/.(t)|P < ^ - yh{t) < j - a{t)Ea/ 2 {b{t)t°^/‘^), 

1 — a 1 — a 

and the result follows after using (3.5) to estimate a{t), because the lower bound 
r(Q!/2 + 1) > 4/5 for 1/2 < a < 1 implies b(t) < |CFt“^^/(l — a). Note that (3.5) 
does not rely on the first assumption Tj > 0 of Theorem 3.2. □ 

To estimate the error in the finite element solution, we will compare Uh (t) to the 
Ritz projection of u{t). Recall that Rh ■ Ho{^) —t is defined by 

{{Rhw)x,Vx) = {wx,Vx) for all w e 
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and satisfies 


\\v - Rhv\\<ChJ'\v\r and \\(v - Rhv)x\\ < ChJ' (3.6) 

for r S {1,2} (in our piecewise-linear case). Here, |u|r = ||u^’’)|| is the usual H '"- 
seminorm. The next theorem shows that if we choose uoh = RhUo, and if u S 
((0,T),i?''(H)), then \\ufi{t) — u{t)\\ = 0{h'~) for 0 < t < T and r £ {1,2}. 
Theorem 3.4. Let Uh denote the the spatially-discrete finite element solution 
of (1.1), defined by (3.2) with gt,{t) = 0. Then, under the hypotheses of Theorem 3.3, 
we have the error bound 

||u/.(t) - u{t)f < C\\uoh - RhUoWl + Ch^^ (^lluoll? + ^ lkt(s)||? ds 

for 0 <t <T and r £ {1,2}, where C depends on a, F, T and L. 

Proof. We decompose the error into two terms, 

Uh — u = 9 p where 9 = Uh — RhU and p = RhU — u, 

and deduce from (3.2) that, for v £Sh, 

{9t,v) + {dl~°‘9,,,va;) - {Fd}~°‘9,v,f) = {g{t),v) 

- {RhUt,v) - {dl~°‘{RhU)a:,Va:) + {Fdl~°‘RhU,Vx). 

Since {dl~°'{Rhu)x,Vj;) = {{Rhd]~°‘u)x,Vx) = {d]~°‘ux,Vx), it follows from (3.1) 
that 0 : [0, T] —>■ Sft, satisfies 

{9t,v) + {dl-^9,,v^) - {Fdt^9,v,) = {Fdl-‘^p,v,) - {pt,v), 

which has the same form as (3.2), with 9, —pt and Fdl~°‘p playing the roles of Uh, 
g{t) and g*{t), respectively. Hence, Theorem 3.3 gives 

\\9{T) - 0(O)f < C||0(O)||{ + C Tdlptf + Wdl-^pf) dt, 

^0 

and by Lemma 2.4, pW^ dt < (711/9(0)11^ + C ||pt|pdt. The desired error 

bound follows after applying (3.6) with v = Ut. □ 

4. An implicit time-stepping scheme. To discretize in time, we suppose 
0 = to < ti < t 2 < ■ ■ ■ < t]s[ = T and denote by kn = tn — t„_i the length of the 
nth subinterval In = (tn-i,tn), for 1 < n < N. The maximum time step is denoted 
by fc = maxi<„</v kn. With any sequence of values v^, ..., we associate the 

piecewise-constant functions v and v defined by 

v{t) = u" and v{t) = for t„_i <t<tn. (4.1) 

Integrating the fractional Fokker-Planck equation (1.1) over the nth time interval /„ 
gives 


i{tn) - u{tn-i) - J dl °‘u„j;dt + j [Fd] °‘u)^dt = j g(t)dt. 


(4.2) 
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We seek to compute U'^{x) « u{x, tn) for n = 1, 2, ..., iV by requiring that 

_ c/"-i _ f dt+ [ dt = Kr, (4.3) 

■If -If 

with F"‘{x) = F{x,tn) and g" Ri h 9 {'t)dt. The time stepping starts from the 
initial condition 


U°{x) = uo(x) for 0 < X < L, 

and is subject to the boundary conditions C/”(0) = 0 = U'^{L) for 1 < n < N. 
Since 


(4.4) 




LUniV-' 


j=l 


i=i 


where 


-'nj 


/ UJa{tn - S) ds = LOi+aiin - tj-l) - lUi+a{tn - tj) for n > 2, 

■hi 


with wii = a;i+Q.(ti), we see that 



d] °^vdt = {I°‘v){tn) - {I°'v){tn-l) = ^UJnjV^ 

i=i 


n—1 

i=i 


Hence, to find C/" satisfying (4.3) we solve 


n—1 

- OJrinK, + a;„„(F"C/"), = 

i=i 


ix 




It follows from Theorem 4.3 below that this linear elliptic boundary-value problem 
has a unique solution 17” € Fl^ih) if k is sufficiently small. Note that if the mesh is 
uniform, that is, if fc = for all n, then the sums are discrete convolutions because 


UJnn — 


Jc^ n 
^ ^n—j 

r(l-ha) 


= uja+i{k)a„-j where On = (n-I-1)“ — : 


(4.5) 


The next two lemmas, which will help prove a stability estimate for C/”, use the 
following notation for the backward difference. 


dv{t) = dv^ 


for t G In- 


Lemma 4.1. For any sequence (v")))Lq in L 2 {Il), 

N N n N 

Y,kn\\{I^dv){tnW<2u^oi+l{T)Y,U^NnY,k,\\dhf+2Y,kl^+^\\dv-^\;^ 

n—1 n—1 j—^ n—1 
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Proof. For t & In, 


||(/“5z;)(i„)|| < / Wait - s)\\dv{s)\\ ds + a;a(tn - s)||9z;(s)|| ds 


= {I°‘\\dv\\){t) + UJa+l{tn - t)\\dv^\\, 

where we used the fact that iOaifn — s) < toa{t — s) because t < tn- Thus, after 
squaring and integrating over t G In, obtain 

kn\\irdv)itn)f = \\irdv)itn)fdt<2j^ ii‘^\\dv\\fdt + 2kl^+^\\dv^f, 

since (2a + l)r(a + 1)^ > 1. By Lemma 2.3, 

N pt 

^ fc„||(/“5w)(t„)|p < 2wa+i(r) / Ua{T-t) \\dv{s)\\'^ ds dt 

-^0 Jo 

+ 2f]e+'l|9' 


,n II2 


n—1 


and the result follows because ||3r;(s)|p ds = kjWdv^W^. □ 

Lemma 4.2. For uniform time steps kn = k and for any sequence (vn)n= 0 ! 


J dl °‘v dt = k(I'^dv)(tn) + ijJnlV° 


and 


N 


N 


n—1 n—1 

Proof It follows from (4.5) that Un-ij-i = ^nj- Thus, 


n — 1 


I°‘v{tn-l) = ^ LOn-ljV^ = ^ Uln-l,j-lV^ ^ ^ UJnjV^ ^ = I°‘v(tn) “ Wnlf° 

j=l j=2 j=2 

and so 


J dl °‘vdt = I°‘v(tn) - I^^vitn-l) = I'^iv - v)(tn) +UJnlV°, 

which gives the first result because v — v = kdv. To prove the second result, use (4.5) 
to write {I°‘v)(tn) = U!a+i{k) cin-jV^ and apply Lemma A.2 (from Appendix A) 

to deduce that, pointwise in x, 

N N n N 

^ V^I^v(tn) = U;a+l(fc) ^ ^ an-jV^V^ > ic.a+l(fc) ^(^;”)^ 

n—1 n—1 j=l n—1 

The desired inequality follows after integrating over 17. □ 
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We are now able to show the following stability estimate. 

Theorem 4.3. Assume 1/2 < a < 1 and consider the implicit scheme (4.3) in 
the case of uniform time steps kn = k. If G i/g (17) fl and 

1 + F{x,tnf‘ + Fx{x,tnf' < Cp for 0 < X < L and I < n < N , 

and if k is sufficiently small, then for 1 < n < N , 


\\U^ -U^r <tnEa 


Cif 


2a 


r(a+ 1) 


■ 1=1 ^ 


t: 


2a-l 


2a- 1 


where Ci = 24Ci;’(l + 2Cf), C 2 = 6(1 + 2Cf) and C3 = (>Cp{\ + 4:Cp). 

Proof Put IT" = f/" — [/°. Since the mesh is uniform and = 0, Lemma 4.2 
implies that 


+ U;o.{t)dt = k{rdW:,:,){tn)+U;nlUl 


XX 1 


and similarly 


^ (T"a/-“i/)^d7 = fc(T"J“(5i7)(7„))^+w„i(F"C7°),. 

Thus, putting $" = 17° — F'^U^, our time-stepping scheme (4.3) implies that 


fcsiT" - fc(/“aw',,)(t„) = [/" - [/"-' - j dt + LOniUi 


L 


7° 

XX 


= kg"- I (f”/“67)^dt + cu„iC7° 
= kr - fc(F"(J“5W)(t„))^ + 


(4.6) 


We take the inner product of both sides of (4.6) with {I°'dW){tn), then integrate 
by parts with respect to x and use the fact that 1T"(0) = 0 = 1T"(T) to arrive at 

fc(5iT", (/“aw)(t„)) + k\\{rdw,){tr.)f 

= k{r, (/“aw)(t„)) + (fcF"(/“5w)(t„) - (/“aw,)(t„)) 

<\k\\rr + \k\\{i<^dw){t^)f 

+ \k-\k\\F^{I^dW){tn)\\+UJnlW\\f + ifc||(/“5M4)(tn)f. 

Since ||4>"|p<(l + ||T"||i^(^))||C/°||f, 

fc(9w", (/“aw)(t„)) + ifc||(/“aiT,)(t„)f < \k\\rf 

+ (1 + wnil^m) {k\\{I‘^dW){tr.)f + k-^u;l,\\UYi), 

SO after summing over n and applying the second part of Lemma 4.2, we see that 


N 


N 


N 


n—1 


n—1 


n—1 


Y, kWYdwYtnY < E + ^Cp E fc||(/“5w)(tjf 

N 

+ 2Cp\\UYlY^~'^nl- 


(4.7) 
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Recall from (4.5) that a;„i = a;o,+i(fc)a„_i with a„ = (n+ 1)“ — and observe that 
a„ < an°‘~^ for n > 1. Thus, for k sufficiently small, 


N 


^ = 


'2a-l -^-1 


n—1 


r(a + l)2 

< 1 




k‘^ 


n—0 
f^ 2 a-l nN-l 


r(a + l)2 


N -1 


1 + (2“-1)2+^(c 


n—2 


T{ay 


■dy<l + 


t 


2a-l 

N -1 

2a-V 


(4.8) 


where, in the final step, we used the assumption that 1/2 < a < 1 and the fact that 

r(a) > 1. 

In a similar fashion, we next take the inner product of (4.6) with dW^ to obtain 


= k{g^,dw^) - {k{F^{I‘^dW)it^))^,dW^) + 

< |A:||rf + ifc||F/‘(/“91T)(t„) + F^{rdW,){U)f 
+ |fc-^a;2j|$S|p + (i + i + i)fc||51T"f 


and hence 


ifc||aiR"f+ fc((/“5w,)(t„),aiT,”) < Ik\\g^^ + Ik-w„,\\‘^>:y 

+ 3fc||F/^|lL(n)||(/“^^M^)(in)f+ 3fc||F-||L||(/“aiT,)(tJf. 

Since ||'i>"p = — F^U^ — F'^U^iy < CFllCl^Hi, after summing over n it follows 

from the second part of Lemma 4.2 that 

N N N 

Y^^j2 k\\dw^^ < 3 ^ k\\rr + E 

n—1 n—1 n—1 

N N 

+ 6Cf J2 k\\{I^dW){t^)f + 6 Cf ^ k\\{I^dW,){t^)\y, 

n—1 n—1 

which, together with (4.7) and (4.8), implies that 

N 

+ iCi ^ k\\{I^dW){t^)\\y 

n—1 

where 


N / ,2a-l \ 

an = c 2Y. k\\r\? +C3iit/°iid 1 + ). 

Hence, by Lemma 4.1, 

y N N 

Y^ < ^An + iCi (w„+i(T) ^ fc2“+i||5iy"f 

n—1 n—1 

N -1 

< ^An + iCia;«+i(T) ^ ujNnY^ + iCi(a;„+i(T)tcjv^ + k^‘^)Y^. 

n—1 
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For k sufficiently small, the term in on the right-hand side is bounded by \Y^ . 
Therefore, because ui^n < k°‘{N — n)““^/r(a). 


D 'LOL ^ 

< An + ~ where Bn = CiUJa+iitN), 


and so 


n—1 

F" < An + BNk ^ UJaitn “ tj)y^ for 0 < 71 < TV. 

1=0 

Thus, by Lemma 2.7, < ANEa[BNt^) = ANEa[Cit^^ /T{a + 1)). Finally, 

• n 

Y,k\\dw^r 
■1=1 

and the result follows. □ 

We can now prove the following error bound, which implies 

||t/'^-77(t„)|| =0(fc“), 

if u is sufficiently regular and if \\g^ — g{t)\\ < Ck°‘ for t G Ij] recall that |7;|r = ||77^’'^||. 

Theorem 4.4. Assume 1/2 < a < 1 and consider the implicit scheme (4.3) 
in the case of uniform time steps kn = k. If F € Loo ((0, T), W/;,(n)), arid if k is 
sufficiently small, then for 0 < tn < T, 



||lT"|p = 


1=1 




1=1 


IT' p plz 

\\U--u{tn)f <cy^ \\gl-g{t)fdt + ce‘^-^ t\ut\ 
.-1 Jin JO 


2 dt 


1=1 


+ Ck^^ r \ut\ldt + Ck^\\uo\\l+Ck^" r WutWldt, 
J k Jo 


where C depends on a, F and T. 

Proof Denote the error at the nth time level by e" = D” — u(tn). Subtracting 
(4.2) from (4.3) yields 


+ ^ (F"a/-“e)^ dt = fcp”, 

where = Pi + P 2 + Ps for 

Pi = 9 "-T [ 9{t) dt, P2= \ [ 5/-“(n - uU dt, 

JIn JIn 

dt. 

Applying Theorem 4.3, with e" and p" playing the roles of D" and p", and noting 
that e° = 0 by (4.4), we see that 


N 

\\e^f<CY,k\\p"f forl<n<iV. 

n—1 


(4.9) 
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Since p'^ = k ^ Jj (g^ — g) dt, we have 


N 


N 


I ^aidn ^n — 1 5 tyj, 

1 UJa{tn - S) - UJa{tn-l - s), 0 < S < tn-1, 


wr-grdt, (4.10) 

n=l n=ldln 

and if we put 

An(s) = • 

and Snjit) = {t — ^ A„(s) ds for t G Ij, then 

kp 2 = I°‘{u-u)a;xitn) - I°‘{u-u)a;xitn-l) = / An{s){u - u)xx{s) ds 

Jo 

71 p pfj n p 

= E/ ^ri{s) Uxxt{t)dtds = '^ 5nj{t)Uxxt{t){t dt. 

j = l dIj J s J Ij 


Hence, 


N ^ N n p p 

E^II^2lP < tEX! / hxxt{t)\\^{t-tj-i)dt 6nj{tfdt 

^ -1 -1»_1^/t 


n=l 

1^ f 


f/ 


(4.11) 


^E / (^-^i-i)l|w^xi(0ll dnjit) dt. 

i=lJIj n =7 A 


We find that 


dnnid) — f ^ai^n ds — 


'„2a-l 


- {tn - t) 


2a-l 


r(a)2(2a- 1) 
and, since 0 < Ula{tn — s) < UJa{tn-l — s) for S < tn-l, 

fc2a-l _ _ ^)2a-l 


for t G In 


dn,n—l{i') ^ / ^a(^n—1 5) ds — 


r(a)2(2a- 1) 


for t e /„_i, 


whereas if 1 < j < n — 2, then the Mean Value Theorem implies that 
'■* ^ (l-a)2 


dnj(t)^ < [ - S)k]'^ ds < 

Jtj-i 


(n - 1 - for tGlj, 


T{ay 


so 


Snj{ty dt < Ck'^°‘ X 


(n-l-j)-"“-V l<J<n-2, 

1, n — 1 < j < n. 


Thus, 


N 

u 


f Snj{tydt<Ck^‘^f2+ E {n-l-j)~^°‘~A<Ck^°‘, 

A ^ n=j+2 d 
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and therefore by (4.11), 


N 

E 

n—1 


N 

n=l * 


(t ^n—l) II dt. 


It remains to deal with + p'^ 2 : where 

P^i = l I -u)+ F-dp^iu - u),) dt, 

P32 = ll {(P - Pn.dl-^^u + {F- F^)dp^u,) dt. 
Estimating p^^ in the same way as p 2 ^ we see that 


N 




sill — 


< Ck^ 


n—1 


p ptN 

/ {t - tn-l)\\ut\\ldt < 0^°" / 

n=l-^F Jo 


WutWfdt. 


(4.12) 


(4.13) 


Next, since \\F{t) — F"||i < Ck for t G In, 

\\P32r<k-^l\\P{t)-PXdtl Wdl-^^uWldt^CkJ^ Wdl-^^uWldt, 

SO, using Lemma 2.4, 

^ ptN / ptN \ 

I \\dl-^u\\ldt<Cei^\\uo\\l+WutWldtj. (4.14) 


The error bound now follows from (4.9), (4.10) and (4.12)-(4.14). □ 

5. Numerical experiments. Our discrete-time solution U" G Hq{I1) of (4.3) 
satisfies 

{dl-^U,,vPdt- {F^dl-^U,vPdt = {g,v)dt 

for all V G Hq{I1). We therefore seek a fully-discrete solution UJ^ G Sh given by 

{u;: -Uj^-\v)+ {dl-^Uh.,vPdt- {F^dl-^Uh,vPdt = {g,v}dt 

for all V £ Sh and for 1 < n < N, with 17° = RhUQ. (In our case, the Ritz projec¬ 
tion RhUo is simply the nodal interpolant to uq.) Explicitly, let (j)p G §/t denote the 
pth nodal basis function, so that 4>p{xq) = Spq and 

p-i 

Uhix) = E where U; = UJi{xp) « U^{xp) « u{xp, t„). 

p=i 

Define the (R—1) x (P—1) tridiagonal matrices M and J5" with entries Mpq = (</>,, ^p) 
and = {(j)qx ; 4^px ) - {F"4lq,(l)p,^), and define (P — l)-dimensional column vectors 
Lf" and G" with components Up and G” = {g, 4>p) dt. We find that 


n n—1 

MG" - MG""^ -f E - E w„-i jB"G^' = G", 


1=1 


1=1 
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Table 5.1 

Behaviour of h, defined by (5.1), as the number of time steps N increases, for different 
choices of the mesh qradinq parameter'y. In each case, a = 0.625 and we use a spatial resolution h = 
L/P with P = 5120. 


N 

7 = 1.0 

rt 

7 = a ^ = 1.6 

rt 

II 

to 

0 

rt 

80 

8.93e-03 


8.60e-03 


l.Ole-02 


160 

4.95e-03 

0.851 

4.33e-03 

0.989 

5.09e-03 

0.986 

320 

2.80e-03 

0.823 

2.18e-03 

0.993 

2.56e-03 

0.992 

640 

1.62e-03 

0.791 

1.09e-03 

0.996 

1.28e-03 

0.995 


Table 5.2 

Behaviour of h, defined by (5.1), as the number of time steps N increases, for different 
choices of a. In each case, 7 = 1 and we use a spatial resolution h = L/P with P = 5120. 


N 

a = 0.25 

rt 

a = 0.50 

rt 

a = 0.75 

rt 

80 

1.93e-01 


2 . 21 e -02 


7.40e-03 


160 

1.70e-01 

0.183 

1.50e-02 

0.554 

3.73e-03 

0.989 

320 

1.50e-01 

0.188 

1.04e-02 

0.538 

1.88e-03 

0.990 

640 

1.31e-01 

0.193 

7.20e-03 

0.525 

9.46e-04 

0.989 


SO at the nth time step we must solve the linear system 

n —1 

(M + 0JnuB^)U^ = MU^-^ + G” - y] 

We now describe some experiments using this numerical scheme. 

5.1. Convergence behaviour. In our first test problem, we considered (1.1) 
with 


F(x,t) = X + sint, T = 1, L = tt, Kq =/Xq, = 1, 

where the source term g was chosen so that u{x, t) = [l+wi+Q,(t)] sinx. It follows that 
Ut = as t —)• O’*', and this singular behaviour is known to be typical [12] for 

the fractional diffusion equation (that is, when the lower-order term in F is absent). 
We employed a uniform spatial grid with h = tt/P, but allowed a nonuniform spacing 
in time by putting 


tn = {n/NyT, where 7 > 1 . 

Thus, 7=1 gives a uniform mesh with k = T/N, but if 7 > 1 then the time 
step is initially ki = T/N^ = 0{k'^) and increases monotonically up to a maximum 
oi k = kN ^ jT/N. Such meshes [13] are commonly used to compensate for singular 
behaviour in the derivatives of u at t = 0. As a measure of the error in the numerical 
solution, we computed 


EN,h = Wh - u{tn)\\L2{Q), (5.1) 

(where the spatial L 2 -iiorm was evaluated via Gauss quadrature) and sought to esti¬ 
mate the convergence rates rt and such that 


EN,h « GiF‘ + C2h 














18 


Kim Ngan Le, William McLean and Kassem Mustapha 


Table 5.3 

Behaviour of defined by (5.1), as the spatial resolution h = LfP increases, for different 

choices of a. In each case, 7 = a~^ and we use N = 10,000 time steps. 


p 

a = 0.25 

Tx 

a = 0.50 

Tx 

a = 0.75 

Tx 

4 

8.43e-02 


8.22e-02 


7.74e-02 


8 

2.97e-02 

1.505 

2.92e-02 

1.495 

2.77e-02 

1.483 

16 

6.21e-03 

2.258 

6.07e-03 

2.264 

5.75e-03 

2.268 

32 

1.50e-03 

2.052 

1.46e-03 

2.054 

1.39e-03 

2.046 

64 

3.47e-04 

2.108 

3.23e-04 

2.177 

3.03e-04 

2.201 


Fig. 5.1. Estimated convergence rate rt as a function of a, with uniform time steps. 








r 







( 








A 

i 







A 

r 

< 






A 

r 

r' 









0 . 0 1-^^^^. . .^-1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 


a 


from the relations 

rt « rt(iV, h) = \og2{EM,h/E2N,h) when -C F*, 

Tx ~ rx{N, h) = \og2{EM,2h/EN,h) when fc’'* <C . 

We first tested the convergence behaviour with respect to the time discretization. 
Table 5.1 shows how Ejsr^h varies with N, for a fixed, high-resolution spatial grid with 
P = 5120 subintervals, when a = 0.625 and for three choices of 7. In the case of a 
uniform mesh (7 = 1), we observe rt « 0.8, suggesting that the 0{k°‘) error bound of 
Theorem 4.4 is somewhat pessimistic in this case. Although the convergence analysis 
of our time-stepping scheme applies only when 7 = 1, we observe that ~ Ck if 
7 > a~^ = 1.6. (The constant C is smallest when 7 = a~^.) Table 5.2 shows results 
for three different choices of a as we vary N, using uniform time steps (7 = 1) and 
the same fixed spatial grid as before. Note that the choices a = 0.25 and a = 0.5 
are excluded by our theory, which requires 1/2 < a < 1. Figure 5.1 gives a more 
complete picture of the convergence rate rt as a function of a when 7 = 1, and may 
be compared with the known result rt = min(2a, 1) for the homogeneous diffusion 
equation (that is, the special case F = 0 and g = 0) with regular initial data [14, 
Lemma 6]. 

Next, we tested how Ejq^h behaves as the spatial mesh is refined, using a fixed, 
high-resolution time discretization with N = 10, 000. Table 5.3 shows results for 
three different choices of a using a mesh grading 7 = 0;“^ in each case. We see that 
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Fig. 5.2. Contour plot of the solution for the problem of Section 5.2. The dashed line shows 
the first moment x{t). 



Fig. 5.3. Top: first moment (as computed via Laplace transforms) of the solution for the 
problem of Section 5.2. Bottom: error in the first moment of UJ(. 




t 


E{N, h) Ri Cih^, consistent with Theorem 3.4 (when 1/2 < a < 1). 

5.2. An application. In our second example, we solve the homogeneous equa¬ 
tion on the spatial interval (—L,i), that is, 

Ut — dl~°‘uxx + X ~ ® for 0 < t < r and —L < x < L, 

with F = —X -\- sint, subject to the boundary conditions m(±L, t) = 0. For the initial 
data uq, we chose a normal probability density function with mean 0 and variance . 
This choice of F is taken from a recent paper by Angstmann et al. [1]; notice that 
Fa; = —1 < 0 so the first assumption of Theorem 3.2 is not satisfied and we must 
rely on Theorem 3.3 to ensure stability of the spatially discrete scheme (3.2). For our 
computations, we used the values a = 0.75 and cr = 0.5, with a mesh grading param¬ 
eter 7 = l/a. Figure 5.2 shows a contour plot of the numerical solution computed 
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using our fully discrete method in the case L = 9, T = 10, N = 100 and P = 2L^. 
Although we do not know an analytical solution, Laplace transform techniques [1] 
show that in the limiting case when L —>■ oo, and interpreting u{-,t) as a probability 
density function, the expected position, or first moment, is 


/ oo If* 

xu{x,t)dx =—,—o / Ea(—(t — s))s°‘~^smsds, 

-oo r(a) Jo 


where Ea denotes the Mittag-Lefller function (2.1). Figure 5.3 shows the oscillatory 
behaviour of x{t) for 0 < t < T = 70, and the difference between this theoretical value 
and the first moment of the numerical solution in the case L = 20, N = 20T and 
P = 2TL^. We observe little if any loss of accuracy over more than 10 oscillations. 

5.3. Non-smooth initial data. In the special case of a fractional diffusion 
equation {P = 0 and g = 0), a standard energy argument shows that both the exact 
solution and the finite element solution are stable in L 2 (n), with 


||u(t)|| < llitoll and ||u?>(t)|| < ||mo?>|| for t > 0 . 


By comparison, for nonzero F the stability estimates of Theorems 3.2 and 3.3 yield 
weaker bounds of the form 


\\uh{t)\\ < C\\uoh\\i for 0 < t < r, (5.2) 

and in the case of our (spatially continuous) time-stepping scheme. Theorem 4.3, 

\\U^\\<C\\U °\\2 for 0 <t<T. (5.3) 

To investigate whether the stability properties really depend on the smoothness of the 
initial data, we solved the same problem as in Section 5.2 but chose the nodal values 
for the discrete initial data to be uniformly distributed pseudorandom numbers in 
the unit interval. For 0 < t < T = 40 and many different combinations of N and P, 
we never observed any kind of instability. In all cases, the solution quickly smoothed 
and began an oscillatory behaviour similar to that seen in Figure 5.2, suggesting that 
(5.2) and (5.3) are pessimistic with respect to the regularity required of the initial 
data. 

Appendix A. Positivity of discrete convolution operators. 

Recall the following positivity property of Fourier cosine series. 

Lemma A.l. If the sequence uq, ai, 02 , ... tends to zero and satisfies 

a„ > 0 and a„+i < 5 ( 0 ^ -|- an+ 2 ) for all n > 0, 


then 

OO 

— + On COS nO >0 for —tt < 6* < tt. 

Proof Zygmund [20, p. 93 and Theorem (1.5), p. 183]. □ 

For 0 < a < 1, let 

n 

{AU)^ = ^ On-jU^ where On = {n + 1 )“ — n“. 
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We used the following inequality in the proof of Lemma 4.2. 

Lemma A.2. For any real, squart-summable sequence , ..., 


OO 


OO 


j2iAuru->-j2{un 

n—0 n—0 

Proof. Define V{0) = and observe that 


OO OO 


u{6)v{e) de = EE 

ri=0 j=0 


J{n-j)e 


de = 


n—0 


Since 


OO / n 


OO / OO 


AU{e) = ^ f ^ j = E ( E = d(0)U(0) 

n—0 '^j—0 ^ j—0 ^n—j 


we conclude 


OO -| pTz ^ _ 

y"{AU)^W=— d{0)U{9)V{0)d9. 

t'o ^TT J_^ 

In particular, when 14” = [/” is purely real, 

OO «7r ^ 

J2{Auru- = - m{9)\U{9)\ 

n=0 


d0. 


The function f{x) = (cc -I- 1)“ — a;“ is positive for a; > 0, and as a: ^ oo, 

f{x) = x“[(l + x-^y - 1] = a;“[aa:-i + 0{x-^)] = ax^-^ + 0(a:“-2), 

so in particular f{x) —> 0. Furthermore, / is convex because 

f"{x) = a{a — 1) [(a; -I- 1)““^ — a;““^] = q;(1 — a) [a;““^ — (a; -I- 1)“”^] > 0 

for all a; > 0, so the sequence a„ = f{n) satisfies the assumptions of Lemma A.l. 
Hence, 3?a(0) > ao/2 = 1/2 and finally 


OO -I pTT -| _ OO 

Y.iAUTU->- -\u{0)\^d0 = -Y,{U 

n=0 ^ n=0 
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